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Abstract. 

In this work, charm production in the DPMJET hadronic jet simulation is 
compared to experimental data. Since the major application of DPMJET is the 
simulation of cosmic ray-induced air showers, the version of the code integrated 
in the CORSIKA simulation package has been used for the comparison. Wherever 
necessary, adjustments have been made to improve agreement between simulation 
and data. With the availability of new muon/neutrino detectors that combine 
a large fiducial volume with large amounts of shielding, investigation of prompt 
muons and neutrinos from cosmic ray interactions will be feasible for the first time. 
Furthermore, above > 100 TeV charmed particle decay becomes the dominant 
background for diffuse extraterrestrial neutrino flux searches. A reliable method to 
simulate charm production in high-energy proton-nucleon interactions is therefore 
required. 



1. Introduction 

During recent years neutrino detectors, notably IceCube and ANTARES, as well as 
large air showers arrays, such as the Pierre Auger Observatory and the Telescope 
Array, have been advancing or completing their construction. These large detectors 
are sensitive to an energy region in which contributions from prompt charm decays in 
cosmic ray showers cannot be neglected and may constitute an interesting signal as well 
as a significant background. In the case of searches for diffuse fluxes of astrophysical 
neutrinos in neutrino telescopes, the signal must to be separated at high energies from 
the background of atmospheric neutrinos. 

Atmospheric neutrinos from 7r and K meson decays, whose spectrum 
approximately follows the power law E~ 3 7 above lOOGeV, can be rejected using cuts 
on energy-correlated variables. However, the spectrum of prompt neutrinos is expected 
to be harder and correspond to the primary cosmic ray spectrum, since charmed 
mesons tend to decay promptly without interacting in the atmosphere (see e.g. dJ[2])- 
Consequently, this contribution needs to be known and separated from the signal, 
which is generally assumed to follow a power law of w E~ 2 , characteristic of Fermi 
acceleration processes [3] . It is therefore essential to achieve a precise understanding of 
the expected charm-induced lepton flux and angular distribution, in order to identify 
it in the data above 100-lPeV, where the contribution of the softer component from 
light meson decays becomes negligible. For this purpose, the entire air shower needs 
to be simulated. 
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One of the most widely used cosmic ray shower simulation software packages in 
high-energy astrophysics is CORSIKA jjj . It offers the user the choice between several 
high-energy hadronic jet simulations but does not yet officially support propagation 
and decay of charmed particles. In order to remedy this situation, first the viability 
of one of the available models for simulation of charm has to be demonstrated. 

The consistency of charm simulation with experiment was first verified in 1995 
using the standalone version of dpmjet[5]. The best available measurements at the 
time were those from the LEBC hydrogen bubble chamber that took data at CERN 
and Fermilab during the late 1980s using proton beams of 400GeV and 800GeV, 
respectively [6l [7] . Since then, a significant amount of new experiments have taken 
place, most notably in fixed-target runs at Fermilab, but also using various types 
of colliders and most recently the HERA-B detector at DESY [8]. In this work, 
the version of DPMJET integrated in CORSIKA was used, running in first-interaction 
mode. This was done to directly evaluate its applicability to air shower simulation 
and to ensure compatibility with the standalone version. 

2. Charm production in the Dual Parton Model 

2.1. The Dual Parton Model 

There are two mechanisms in any model for hadron production: soft particle 
production and hard collisions. Both components are part of the DPMJET hadron 
production models used in this work. Soft hadron production cannot be rigorously 
derived from QCD, the gauge field theory of strong interactions. Most theoretical 
efforts use Regge theory for a systematic description of soft hadron production. The 
so-called hard component of particle production can be derived from perturbative 
QCD. QCD perturbation theory within the QCD improved parton model can only 
be applied at moderately large transverse momenta (p±). One feature of soft 
hadron production at small p± values is represented by exponentially decreasing 
transverse momentum distributions. The QCD improved parton model predicts 
another component of the transverse momentum distributions, decreasing like a power 
law, and thus less steeply than an exponential function, at large p±. 

The Dual Parton Model DPM [9] and the closely related Quark Gluon Strings 
Model QGSM [TO] ITTj construct hadron production from fragmenting strings. Here 
we will describe the structure of this mechanism in hadron-hadron collisions. Since it 
is not possible to describe soft component of hadron physics using perturbative QCD, 
't Hooft introduced a new expansion parameter [12]. First, QCD is generalized from 
the gauge group SU(3) to SU(iV c ) with N c representing the number of quark "colours" . 
The idea was to use 1/N C as an artificial expansion parameter, and later set N c = 3 for 
physical applications. One then finds that the Feynman graphs can be characterized 
by a two-dimensional surface, and the Feynman expansion can be regrouped as a 
sum over surface topologies. The leading order corresponds to planar graphs, the next 
order for 2-to-2 amplitudes involves a cylindrical topology. It has been shown [13] that 
there is a one-to-one correspondence between the terms of this topological expansion 
and the terms of Reggeon Field Theory RFT[T4]. The second theoretical concept 
entering DPM is duality [15] and the Dual Topological Unitarization DTU[16]. The 
third ingredient of DPM is the coloured parton model. 

The leading contribution to soft hadron production in the DPM corresponds to 
the production of two chains. The two chains arise from the unitarity cut of the 
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Pomeron exchange diagram. Located at the end of either chain are the valence quarks 
and diquarks of the two protons. The valence diquarks have anti-triplet colour, while 
the colliding protons as well as the two produced chains are colour neutral. Particle 
production occurs via fragmentation of the two quark-diquark chains. In order to 
arrive at a quantitative model, we have to specify the probability pi(x) that the 
interaction separates the proton into a valence quark with momentum fraction x and 
a valence diquark with the remaining momentum fraction (1 — x). These longitudinal 
momentum fractions are given by Regge asymptotics. Valence quarks in baryons as 
well as in mesons follow a;^ 1 / 2 , valence diquarks in baryons (1 — a;) 3 / 2 and sea quarks 
x^ 1 distributions. Combining these, we obtain for the distribution in a proton 

p\ = clx~^{l-xf/\ (1) 

and in a meson (valence quarks follow the same distribution as their corresponding 
antiparticles) 

p™ = c ?x- 1/2 (l-x)- 1/2 . (2) 

where the coefficients c\ and c™ normalize the distributions to unity. 

This two-chain model already describes hadron production at low energies quite 
successfully. Many Monte Carlo models therefore started with two-chain models. 
However, at higher energies more and more discrepancies between the model and 
experimental data appear. One example is that the measured transverse momentum 
distributions can only be explained by including a hard component in the hadron 
production. Also, from a theoretical point of view, the pomeron total cross section 
(using the so-called supercritical pomeron) violates the unitarity principle. Therefore, 
the model has to be unitarized. 

The DTUJETpT] and PHOJETpS [19] Monte Carlo models use an eikonal 
unitarization scheme, dtujet and phojet, respectively, are used for the elementary 
interactions in DPMJET-II and dpmjet-III. Both models have multiple soft and hard 
chains as demanded by the unitarization method. 

DPMJET-II [3D] and DPMJET-III [5T] use the Glauber model to describe hadron- 
nucleus and nucleus-nucleus collisions. Here, we will not present any details about the 
Glauber model and refer the reader instead to the original publications cited above. 

2. 2. Charm production in DPMJET and CORSIKA 

The paper by Battistoni et al. [5] gives the first description of DPM charm production. 
They used a version interfaced to the cosmic ray cascade code HEMAS [22], which 
has since then been phased out. There, DPM was implemented in form of the Monte 
Carlo event generator DPMJET-II [5T] , which itself is based on earlier work [531 (Ml US] • 
In [5], hard charm production was tested against NLO perturbative calculations 
|26j and differential charm production cross sections were compared to experimental 
data, mainly from the LEBC-EHS and LEBC-MPS collaborations 0[7]. Hadronic jet 
simulation with dpm jet-II is available in the CORSIKA cosmic ray transport code [4] . 
However, the official CORSIKA version at the time DPMJET-II was implemented had 
no provision for the propagation of charmed hadrons. Instead, all charmed hadrons not 
decaying within the hadronic jet simulation itself were handled by substituting strange 
quarks for charm. Charm production from DPMJET-II and other jet simulations 
is expected to become available in the official CORSIKA release [29] in the near 
future. Earlier, an unofficial modified version was provided to the authors [57], in 
which charm was no longer supressed. Our work is mostly based on this version, in 
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which however one crucial parameter was erroneously changed with respect to the 
original DPMJET. The comparisons with experiment shown below demonstrate both 
the influence of correcting this parameter and the modifications made subsequently 
to the main DPMJET code (see Figs. [3] and 2]). 

dpmjet-II is also integrated in the FLUKA hadron cascade code [28 . We are 
not aware of any charm calculations using FLUKA-dpmjet-II, but this might change 
soon. All modifications to charm production described in this paper are similarly 
applicable to FLUKA. 



2.3. Charm production in elementary hadronic collisions 

In dpmjet-II there are three different mechanisms by which production of heavy 
flavors like charm can occur: (i) charmed quarks produced at the ends of hard and 
semihard chains (minijets), (ii) charmed quarks produced at the ends of soft sea chains 
and (iii) charm production inside the chain decay. Only the first mechanism (i) is 
founded on the solid theoretical basis of perturbative QCD, the other two mechanisms 
are phcnomcnological models used for charm and prompt muon production within 
cosmic ray-induced cascades 



2.4- Charm production at the ends of hard and semihard chains (minijets) 

Regarding this component, nothing in dpmjet-II was changed with respect to the 
older version described in [5]. Therefore, we do not need to repeat the details given 
there. However, it should be mentioned that in that work a detailed comparison 
was presented between the results of the dpmjet-II Monte Carlo at large transverse 
momenta and the corresponding NLO QCD calculations, showing excellent agreement 
between the two. For hadronic and nuclear collisions at laboratory energies beyond 
the TeV range, this is the dominant mechanism of charm production. 



2.5. Charm production at the ends of soft sea chains 

In DPMJET, charm production at the end of soft sea chains is considered as the non- 
perturbative limit of minijets. The parton transverse momentum distributions of 
minijets and soft sea chains are joined smoothly at the threshold transverse momentum 
between the two. A certain fraction of the soft sea chain ends carry heavy flavors. 
At high energies, the probability P cS for production of a cc sea quark pair approaches 
the corresponding value for semi-hard chains, while at low energies P cS decreases as 
required by phenomenological considerations, dpmjet uses 

P cS = C / 2Exe- b " E± dEx (3) 

J m q 

with 

&c = 6+1.3 -log 10(|g^). (4) 

An additional constraint is given by the requirement that P cS should not exceed the 
corresponding probability for minijets. 
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2.6. Charm production inside the soft chain fragmentation 

DPMJET-II uses the Lund model PYTHIA for chain fragmentation 30J. Within this 
model, a quark- antiquark pair leading to string breakup is produced via quantum 
mechanical tunneling. The tunneling probability can be worked out as function of 
the transverse mass mj_ of the pair [30j . with the result that cc pair production is 
highly supressed: uu:dd:ss:cc =1:1: 0.3 : 1CP 11 . There are other models for 
string fragmentation with larger cc production probabilities, but none in which charm 
production inside the soft chain could not be neglected. 

In DPMJET-II the picture remains the same, with one exception: the cc probability 
near a valence diquark at the end of a diquark-quark or diquark-antidiquark chain 
is not negligible. This effect was demonstated in hadroproduction by fixed target 
experiments such as SELEX [3U |32]. The leading particle shares a valence quark 
with the incident hadron, while the non-leading one does not. Therefore leading 
particles are copiously produced at large ij? in the forward region of the incident 
hadron, resulting in an asymmetry between leading and non-leading particles. 
This phenomenon is referred to as leading quark effect [33j . SELEX finds in the 
fragmentation region of protons and E _ hyperons A+ charmed hyperons with a 
comparatively flat Feynman-x (xf) distribution. In a fit to (1 — |:ef|) q , ol is about 2.5. 
Without a new mechanism, dpmjet generates this distribution with an a of about 7. 
Better agreement with SELEX can be achieved by introducing, adjacent to a diquark, 
cc pairs with increased probability. In order to retain the standard pythia code, this 
change was done in a special routine of dpmjet-II. The comparison of SELEX data 
and the modified version of dpmjet is shown in section 3.3. 



3. Comparison with Experimental Data 

3.1. Collision Scaling of Charm Production 

In dpmjet, nuclear collisions are treated with a Monte-Carlo formulation of the 
Glauber model [201 121j. This model contains two important numbers: the number 
of participants N part and the total number of collisions N co u. Soft particle production 
scales rather well with N part , whereas hard particle production (for instance at large 
p±) scales with N co u. The latter phenomenon is commonly referred to as collision 
scaling. All charm production, because of the large mass difference of the quarks 
involved, is expected to follow collision scaling. 

This behaviour can be reproduced in experimental data. Fermilab experiment 
E789 [33] measured neutral D-meson production in p-Be and p-Au collisions, finding 
the production cross section to behave like 

<t(A) = a A a (5) 

with a = 1.02 ± -.03 ± 0.02 and the ratio 

R= a{ j^J = 1.06 ±0.11 ±0.07. (6) 

The consistency of these values with unity corresponds to collision scaling. 

The HERA-B Collaboration [S] measured D°, D+, D+ and D*+ in p-C, p-Ti and 
p-W collisions. They parametrize the production cross sections as 

a p -A — o- p _ N A a (7) 

and obtain a — 0.99 ± 0.04 ± 0.03 in agreement with collision scaling. 
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Table 1. Comparison of collision scaling in dpmjet-II with the results from 
HERA-B [5]. (all cross sections in mbarn) 
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Figure 1. Total open charm production cross section measurements, taken from 
1 351 and references therein, represented by markers. The corresponding values 
from DPMJET simulation are indicated by the dotted line. 



DPMJET-II in its original form included charm production in which collision 
scaling was represented rather poorly. This has been corrected in the version used 
in this paper. 

In Table Q] we compare the o~ p -n cross sections obtained from dpmjet-II in p- 
p and p-Nitrogen collisions with the ones obtained by the HERA-B collaboration. 
Taking into account experimental errors (as given in the table), we find good 
agreement, meaning that collision scaling as simulated in dpmjet-II corresponds to 
that found experimentally. 

3.2. Total Charm Production Cross Section 

Figure [1] shows the total open charm production cross section as measured, along with 
the corresponding values from dpmjet. It can be seen that the simulation reproduces 
the actual values reasonably well, considering that the errors on the experimental 
values are quite large. 



Charm Production in DPMJET 



7 



LEBC-EHS (all D) 



LEBC-EHS (all D) 



> 

a 10 k 

"a 



DPMJET(old) 
— DPMJET(new) 
. — *— Experimental 



+ 



1 2 



LEBC-MPS (all D) 



> in 1 — 

o 
O 

? I 



8 9 ,10 
p 2 [GeV 2 /c 2 ] 



— DPMJET(old) 

— DPM.IET(new) 
— *— Experimental 



5 , 
p 2 [GeV 2 /c 2 l 



10" 



— DPMJET(old) 

— DPMJET(new) 
— *— Experimental 



"i 



T 



i 



i 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

X F 

LEBC-MPS (all D) 



— - DPMJET(old) 
— DPMJET(new) 
— *— Experimental 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 I 



Figure 2. LEBC hydrogen bubble chamber measurements, sum over all D- 
mesons. Values are taken from [B] and [7]. Top left: 400GeV p 2 ± , Top right: 
400GeV x F , Bottom left: 800GeV p 2 ± , Bottom right: 800GeV x F . In this 
case, results from the corrected version (cf. Figs. \3\ and |4j are identical to 
DPMJET(new). 



3.3. Differential Cross Sections 

Precision measurement of differential charm production cross sections at low transverse 
momentum p± and Feynman-x xf = VhlVh mau: requires comprehensive instrumental 
coverage of the interaction region, especially in direction of the beam. Since this 
is difficult to achieve in collider experiments, relevant data are almost exclusively 
provided by fixed-target setups. In the past, various experiements have been conducted 
using different types of particle beams. In the following, we will restrict ourselves to the 
discussion of results that were obtained with proton beams. Even though statistics for 
proton- nucleus interactions are relatively low compared to those for meson beams |37) , 
these are the processes that are the most relevant in cosmic ray-induced air showers. 
Moreover, we consider data in the low p± region most relevant for air showers. 

Figure [5] shows the simplest case, a proton beam impinging on hydrogen nuclei in 
the LEBC bubble chamber. These results were already used in the first comparison 
between DPMJET and experimental data in 1995 5J. Taking into account the limited 
statistics going into the experimental result, agreement between experiment and 
simulation is satisfactory. 

An important aspect of charm hadroproduction is the leading quark effect 
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Figure 3. SELEX measurements of A c production parameters and asymmetry 
for a 540GeV proton beam on a mixed C/Cu target. Top left: Aj!~ xp, Top right: 
xp, Bottom left: p 2 ^ for all A c , Bottom right: Asymmetry between A~i~ and 
A~~ production as a function of xp. Asymmetry values higher than one are the 
result of the likelihood method used to analyze the data. All measurements taken 
from |36| . The error bars have been extended with respect to the original paper 
to reflect uncertainties in absolute A c production cross section, for which a value 
of 7 ± 5.25/^barn/nucleon was assumed (see text). Shown are the result from 
the old release version of CORSIKA, the version with the corrected parameter 
(bugfix) and the newest version in which leading quark effects are implemented. 
For more information about the different versions, see Section 12.21 For the A~f , 
the parameter-corrected result is identical to the new version, in the asymmetry 
plot it is identical to the old version. 



described in section 12.61 by which production of hadrons which share a common 
diquark with the incoming projectile is favored. The extent of this effect has been 
demonstrated by the SELEX colaboration using beams of protons, ir~ and S - . In 
the forward Xf region, production of was found to be significantly enhanced over 
its antiparticle for p and E~ beams, whereas for the ir~ beam no asymmetry was 
observed. Including this effect in the latest vesion of dpmjet leads to good agreement 
with the experimental result, as shown in Figure [3J 

Since the original publication only gave absolute event numbers [5S] , and no value 
for the charm production cross section was published by the collaboration, in Figure 
[3] an additional systematic error was added to the experimental result in order to 
account for uncertainties in the estimate. In 0, the value for A c production was 
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Figure 4. HERA-B measurement for D-meson fixed target production using a 
920GeV proton beam on C, Ti and W targets. Left: xp, Right: p±. Experimental 
data are taken from [B]. For explanation of the different DPMJET results, see 
Figure |3] and Section HOI 



found to be 

lAnbarn < a(pp -> A+/A~)B(A C ^charged) < 6.1^ibarn (8) 
which together with the fraction of A c to total open charm 

, K/K ) = .28 ±0.21 (9) 
cc 

and the total open charm production cross section from [37j 

o~(pp — > cc) = 2b^barn (10) 

results in an estimate for the charm production at SELEX energies (540GeV) as 
7 ± b.25fibarn/nucleon. This is consistent with the value obtained from simulating 
charm production in DPMJET (9.5/ibarn/nucleon). 

Finally, the simulation was compared to data from HERA-B [8]. Even though 
the kinematic region covered by the detector was limited to —0.15 < xj < 0.05 and 
p± < 3.5GeV/c, the result represents the most accurate measurement of open charm 
production using a proton beam. Figure [3] shows the comparison with our simulation. 
Even though there appear to be slight discrepancies in the individual data points, 
the result from DPMJET agrees well with the fit to the data presented by the HERA 
collaboration in their publication in summer 2007. 

SELEX used a mixture of copper and carbon targets, and HERA-B alternated 
between carbon, titanium and tungsten. In each case, the target thickness was a 
small fraction of the proton interaction length. To simplify the simulation, in DPMJET 
nitrogen nuclei were used as target material. This is legitimized by the fact that, 
as noted above, collision scaling was verified both in experiment and simulation. The 
individual diagrams have been scaled to represent the cross section values per nucleon. 

4. Conclusion 

Charm production in dpmjet has been checked against experimental data from 
both fixed target and collider experiments. Since the primary purpose of dpmjet 
is simulation of cosmic ray-induced air showers, the comparison was restricted to 
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proton- nucleoli interactions at low p± , which however represent the by far dominant 
contribution to charm production. The total open charm cross section is consistent 
with experimental values, even though systematic errors on measurements are very 
large, especially at high energies. 

Comparison of simulation to differential distributions from fixed-target 
experiments show reasonable agreement. Since 1995, new data have become available 
that show hard collisions playing a much larger role than previously assumed. Those 
results, most notably from the SELEX collaboration concerning charmed baryon 
asymmetries, but also collision scaling in cc production, were incorporated into the 
code. 

In summary, charm production in DPMJET has been updated reflecting the latest 
experimental results and is now consistent with all available measurements. All 
simulations have been performed with the dpmjet version integrated in the CORSIKA 
air-shower simulation package. It thus allows for the first time to generate Monte-Carlo 
simulations of the prompt component in cosmic ray air showers. The precise effect of 
prompt contributions on muon and neutrino fluxes, as well as a study of meson-nucleon 
collisions, will be subject of a later paper. 
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